i = 0
gbm.results <- list()
ntree_seq <- seq(10,50,500)
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Generate dependent variable
  Violence.mean <-  aggregate(x = as.matrix(dta[dta[,t.var]< year,
                                                dv]),
                              by = list(dta[dta[,t.var]< year,
                                            loc.var]),
                              FUN = mean)
  colnames(Violence.mean) <- c(loc.var,
                               "dv_mean")
  Violence.new <- merge(as.matrix(dta[dta[,t.var]< year,
                                      c(loc.var,t.var,dv)]),
                        Violence.mean,
                        by = loc.var)
  Violence.past <- Violence.new[,dv]-Violence.new$dv_mean
  # Fit GBM
  gbm.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- nrow(dta.past)
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("gbm")) %dopar% {   # Loop over CV runs
                                         set.seed(52*cv) 
                                         shuffle <- sample(1:N,
                                                           N,
                                                           replace=FALSE)
                                         predictions <- matrix(NA,nrow=N,ncol=length(ntree_seq))
                                         sampleSplit = c(0,round((N/5*c(1:5))))
                                         for (f in 1:5) {
                                           # Fit lasso
                                           cv.fit = gbm.fit(x = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                        rhs]),
                                                            y = as.matrix(Violence.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])]]),
                                                            n.trees = max(ntree_seq),
                                                            distribution = gbm_demean_params$distribution,
                                                            interaction.depth = gbm_demean_params$interaction.depth,
                                                            n.minobsinnode = gbm_demean_params$n.minobsinnode,
                                                            verbose = gbm_demean_params$verbose,
                                                            shrinkage = gbm_demean_params$shrinkage,
                                                            train.fraction = gbm_demean_params$train.fraction,
                                                            bag.fraction = gbm_demean_params$bag.fraction)
                                           predictions[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                     ] <- predict(cv.fit,
                                                                  newdata = as.matrix(dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                                    rhs]),
                                                                  type = "response",
                                                                  n.trees = ntree_seq)
                                         }
                                         cv_results <- c("CV Run"=0,
                                                         "Trees"=0,
                                                         "MSE"=1000)
                                         # Get error rates for each lambda used
                                         for (t in 1:length(ntree_seq)) {
                                           if (is.na(predictions[1,t])==0) {
                                             ###!####
                                             MSE <- mean((as.matrix(Violence.past)
                                                        -predictions[,t])^2)
                                             ### ! ####
                                             if (MSE<cv_results["MSE"]) {
                                               cv_results <- c("CV Run"=cv,
                                                               "Trees" = ntree_seq[t],
                                                               "MSE" = MSE)
                                               best.t <- t
                                             }
                                           }
                                         } 
                                         return(cv_results)
                                       }
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  
  gbm.results[[i]]$trees <- round(mean(cv_results[,"Trees"]))
  
  mod = gbm.fit(x = as.matrix(dta[dta[,t.var]< year,
                                       rhs]),
                y = Violence.past,
                n.trees = gbm.results[[i]]$trees,
                distribution = gbm_demean_params$distribution,
                interaction.depth = gbm_demean_params$interaction.depth,
                n.minobsinnode = gbm_demean_params$n.minobsinnode,
                verbose = gbm_demean_params$verbose,
                shrinkage = gbm_demean_params$shrinkage,
                train.fraction = gbm_demean_params$train.fraction,
                bag.fraction = gbm_demean_params$bag.fraction
  )
  gbm.results[[i]]$fit.oos <- predict(mod,
                                      newdata = as.matrix(dta[dta[,t.var] == year,
                                                                   rhs]),
                                      type = "response",
                                      n.trees = gbm.results[[i]]$trees)
  Violence.new <- merge(as.matrix(dta[dta[,t.var]== year,
                                           c(loc.var,t.var,dv)]),
                        Violence.mean,
                        by = loc.var)
  gbm.results[[i]]$actual.oos <- Violence.new[,dv]-Violence.new$dv_mean
  
  print(year)
}
save(gbm.results,file=paste(modeldir,"/",
                              country,
                              "_gbm_",
                              "demean",
                              "_",
                              rhs.group,fileext,
                              ".RData",
                              sep = ""))